package elevation

import (
	"bufio"
	"fmt"
	"image"
	"math"
	"os"
	"path/filepath"
	"strconv"
	"strings"
	"sync"
	"time"

	tiff "github.com/chai2010/tiff"
)

// GeoTransform represents the affine transformation from pixel coordinates to geographic coordinates
type GeoTransform [6]float64

// GeoTIFF represents a GeoTIFF file with its metadata and image data
type GeoTIFF struct {
	Image        image.Image
	GeoTransform GeoTransform
	Width        int
	Height       int
	// 移除缩放因子，因为我们应该直接从数据中获取准确值
}

// openFile opens a file for reading
func openFile(filename string) (*os.File, error) {
	return os.Open(filename)
}

// LoadGeoTIFF loads a GeoTIFF file and extracts its metadata
func LoadGeoTIFF(filename string) (*GeoTIFF, error) {
	// Open the TIFF file
	file, err := openFile(filename)
	if err != nil {
		return nil, fmt.Errorf("failed to open file: %v", err)
	}
	defer file.Close()

	// Get file size for progress estimation
	fileInfo, err := file.Stat()
	if err != nil {
		fmt.Printf("Warning: Could not get file info for progress tracking: %v\n", err)
	} else {
		fmt.Printf("Loading GeoTIFF file: %s (%.2f MB)\n", filename, float64(fileInfo.Size())/(1024*1024))
	}

	// Show loading message
	fmt.Println("Loading GeoTIFF file, please wait...")

	// Start time for performance tracking
	startTime := time.Now()

	// Create a channel to signal when decoding is complete
	done := make(chan bool)
	var img image.Image
	var decodeErr error

	// Show a simple progress indicator since we can't get real progress from the tiff library
	go func() {
		ticker := time.NewTicker(3 * time.Second)
		defer ticker.Stop()
		counter := 0
		for {
			select {
			case <-ticker.C:
				counter++
				fmt.Printf("Loading in progress... (%d seconds)\n", counter*3)
			case <-done:
				return
			}
		}
	}()

	// Decode the TIFF image using the chai2010 library which supports BigTIFF
	// This is run in a goroutine so we can show progress
	var wg sync.WaitGroup
	wg.Add(1)
	go func() {
		defer wg.Done()
		img, decodeErr = tiff.Decode(file)
		done <- true
	}()

	// Wait for decoding to complete
	wg.Wait()

	if decodeErr != nil {
		return nil, fmt.Errorf("failed to decode TIFF: %v", decodeErr)
	}

	// Calculate loading time
	loadingTime := time.Since(startTime)
	fmt.Printf("GeoTIFF loaded successfully in %.2f seconds\n", loadingTime.Seconds())

	// Create GeoTIFF struct
	geoTiff := &GeoTIFF{
		Image:  img,
		Width:  img.Bounds().Dx(),
		Height: img.Bounds().Dy(),
	}

	// Print image dimensions for debugging
	fmt.Printf("Image dimensions: %dx%d\n", geoTiff.Width, geoTiff.Height)
	fmt.Printf("Image type: %T\n", img)

	// Try to load geotransform from TFW file
	geoTiff.GeoTransform, err = LoadTFWFile(filename)
	if err != nil {
		// If TFW file is not available, use default values
		// In a real implementation, these would be parsed from the GeoTIFF tags
		geoTiff.GeoTransform = GeoTransform{
			0, 1, 0, // GT(0), GT(1), GT(2)
			0, 0, -1, // GT(3), GT(4), GT(5)
		}
		fmt.Printf("Warning: Could not load TFW file, using default geotransform values\n")
	} else {
		fmt.Printf("Loaded geotransform from TFW file: %+v\n", geoTiff.GeoTransform)
	}

	return geoTiff, nil
}

// LoadTFWFile loads geotransform parameters from a TFW file
func LoadTFWFile(tifFilename string) (GeoTransform, error) {
	// Replace .tif extension with .tfw
	tfwFilename := strings.TrimSuffix(tifFilename, filepath.Ext(tifFilename)) + ".tfw"

	// Check if TFW file exists
	if _, err := os.Stat(tfwFilename); os.IsNotExist(err) {
		return GeoTransform{}, fmt.Errorf("TFW file not found: %s", tfwFilename)
	}

	// Open TFW file
	file, err := os.Open(tfwFilename)
	if err != nil {
		return GeoTransform{}, fmt.Errorf("failed to open TFW file: %v", err)
	}
	defer file.Close()

	// Read lines from TFW file
	scanner := bufio.NewScanner(file)
	var lines []string
	for scanner.Scan() {
		lines = append(lines, scanner.Text())
	}

	if err := scanner.Err(); err != nil {
		return GeoTransform{}, fmt.Errorf("error reading TFW file: %v", err)
	}

	// Check if we have enough lines
	if len(lines) < 6 {
		return GeoTransform{}, fmt.Errorf("TFW file does not contain enough parameters")
	}

	// Parse the geotransform parameters
	// TFW file format:
	// Line 1: A - pixel size in the x-direction
	// Line 2: D - rotation about y-axis
	// Line 3: B - rotation about x-axis
	// Line 4: E - pixel size in the y-direction
	// Line 5: C - x-coordinate of the center of the upper left pixel
	// Line 6: F - y-coordinate of the center of the upper left pixel

	// Parse values
	A, err := strconv.ParseFloat(strings.TrimSpace(lines[0]), 64) // pixel size in the x-direction
	if err != nil {
		return GeoTransform{}, fmt.Errorf("error parsing TFW parameter A: %v", err)
	}

	D, err := strconv.ParseFloat(strings.TrimSpace(lines[1]), 64) // rotation about y-axis
	if err != nil {
		return GeoTransform{}, fmt.Errorf("error parsing TFW parameter D: %v", err)
	}

	B, err := strconv.ParseFloat(strings.TrimSpace(lines[2]), 64) // rotation about x-axis
	if err != nil {
		return GeoTransform{}, fmt.Errorf("error parsing TFW parameter B: %v", err)
	}

	E, err := strconv.ParseFloat(strings.TrimSpace(lines[3]), 64) // pixel size in the y-direction
	if err != nil {
		return GeoTransform{}, fmt.Errorf("error parsing TFW parameter E: %v", err)
	}

	C, err := strconv.ParseFloat(strings.TrimSpace(lines[4]), 64) // x-coordinate of the center of the upper left pixel
	if err != nil {
		return GeoTransform{}, fmt.Errorf("error parsing TFW parameter C: %v", err)
	}

	F, err := strconv.ParseFloat(strings.TrimSpace(lines[5]), 64) // y-coordinate of the center of the upper left pixel
	if err != nil {
		return GeoTransform{}, fmt.Errorf("error parsing TFW parameter F: %v", err)
	}

	// Map to GeoTransform format according to GDAL convention:
	// GT(0): x-coordinate of the top-left corner of the top-left pixel
	// GT(1): pixel width (x-direction)
	// GT(2): rotation, 0 if image is "north up"
	// GT(3): y-coordinate of the top-left corner of the top-left pixel
	// GT(4): rotation, 0 if image is "north up"
	// GT(5): pixel height (y-direction, typically negative value)
	var gt GeoTransform
	gt[0] = C // x-coordinate of the top-left corner of the top-left pixel
	gt[1] = A // pixel width
	gt[2] = B // rotation term
	gt[3] = F // y-coordinate of the top-left corner of the top-left pixel
	gt[4] = D // rotation term
	gt[5] = E // pixel height

	return gt, nil
}

// LatLonToPixel converts geographic coordinates to pixel coordinates with higher precision
func (g *GeoTIFF) LatLonToPixel(lat, lon float64) (float64, float64, error) {
	// Using inverse of the geotransform equations:
	// X_geo = GT(0) + X_pixel*GT(1) + Y_line*GT(2)
	// Y_geo = GT(3) + X_pixel*GT(4) + Y_line*GT(5)

	// For a north-up image (GT(2) and GT(4) are zero):
	// X_pixel = (X_geo - GT(0)) / GT(1)
	// Y_line = (Y_geo - GT(3)) / GT(5)

	gt := g.GeoTransform
	if gt[2] == 0 && gt[4] == 0 {
		// Simple case - no rotation
		// 使用高精度计算
		xPixel := (lon - gt[0]) / gt[1]
		yLine := (lat - gt[3]) / gt[5]
		return xPixel, yLine, nil
	}

	// General case with rotation - solve the system of equations
	// X_geo = GT(0) + X_pixel*GT(1) + Y_line*GT(2)
	// Y_geo = GT(3) + X_pixel*GT(4) + Y_line*GT(5)

	// Rearranging:
	// X_geo - GT(0) = X_pixel*GT(1) + Y_line*GT(2)
	// Y_geo - GT(3) = X_pixel*GT(4) + Y_line*GT(5)

	// This is a system of linear equations that can be solved using Cramer's rule
	// 使用高精度计算行列式
	det := gt[1]*gt[5] - gt[2]*gt[4]
	if det == 0 {
		return 0, 0, fmt.Errorf("singular transformation matrix")
	}

	xGeoOffset := lon - gt[0]
	yGeoOffset := lat - gt[3]

	// 使用高精度计算
	xPixel := (xGeoOffset*gt[5] - gt[2]*yGeoOffset) / det
	yLine := (gt[1]*yGeoOffset - xGeoOffset*gt[4]) / det

	return xPixel, yLine, nil
}

// GetElevation retrieves the elevation at the specified pixel coordinates using bilinear interpolation
func (g *GeoTIFF) GetElevation(x, y float64) (float64, error) {
	// Check bounds
	if x < 0 || y < 0 || x >= float64(g.Width) || y >= float64(g.Height) {
		return 0, fmt.Errorf("coordinates out of bounds")
	}

	// Use bilinear interpolation
	return g.bilinearInterpolation(x, y)
}

// bilinearInterpolation performs bilinear interpolation to get the value at non-integer coordinates
func (g *GeoTIFF) bilinearInterpolation(x, y float64) (float64, error) {
	// Get the four surrounding pixels
	x0 := int(math.Floor(x))
	y0 := int(math.Floor(y))
	x1 := x0 + 1
	y1 := y0 + 1

	// Check bounds
	if x0 < 0 || y0 < 0 || x1 >= g.Width || y1 >= g.Height {
		return 0, fmt.Errorf("coordinates out of bounds")
	}

	// Get the values at the four corners
	q00, err := g.getPixelValue(x0, y0)
	if err != nil {
		return 0, err
	}
	q01, err := g.getPixelValue(x0, y1)
	if err != nil {
		return 0, err
	}
	q10, err := g.getPixelValue(x1, y0)
	if err != nil {
		return 0, err
	}
	q11, err := g.getPixelValue(x1, y1)
	if err != nil {
		return 0, err
	}

	// Calculate the weights with higher precision
	// 使用math.Floor和math.Ceil确保精度
	dx := x - math.Floor(x)
	dy := y - math.Floor(y)

	// Perform bilinear interpolation with explicit steps for clarity
	// f(x,y) = (1-dx)(1-dy)*q00 + dx(1-dy)*q10 + (1-dx)dy*q01 + dxdy*q11
	// 分步计算以提高精度
	f1 := (1 - dx) * (1 - dy) * q00
	f2 := dx * (1 - dy) * q10
	f3 := (1 - dx) * dy * q01
	f4 := dx * dy * q11
	result := f1 + f2 + f3 + f4

	return result, nil
}

// getPixelValue retrieves the elevation value at the specified pixel coordinates
func (g *GeoTIFF) getPixelValue(x, y int) (float64, error) {
	// For a grayscale image, we can get the value directly
	grayImg, ok := g.Image.(*image.Gray16)
	if ok {
		// For Gray16 images, the value is directly available
		value := grayImg.Gray16At(x, y)
		// 直接返回像素值作为高程值
		return float64(value.Y), nil
	}

	// Try Gray image as fallback
	grayImg8, ok := g.Image.(*image.Gray)
	if ok {
		// For Gray images, the value is directly available
		value := grayImg8.GrayAt(x, y)
		return float64(value.Y), nil
	}

	// For other image types, we'll convert to grayscale
	color := g.Image.At(x, y)
	r, gVal, b, _ := color.RGBA()

	// Convert to grayscale using standard luminance formula
	// This is a simplification - in a real implementation, we would need to handle
	// the actual data type of the DEM (which might be 16-bit or 32-bit float)
	luminance := 0.299*float64(r) + 0.587*float64(gVal) + 0.114*float64(b)

	// Normalize to 0-255 range
	elevation := luminance / 256.0

	return elevation, nil
}
